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We have used diffusion Monte Carlo (DMC) calculations to study the structural properties of 
magnesium hydride (MgH2), including the pressure-volume equation of state, the cohesive energy 
and the enthalpy of formation from magnesium bulk and hydrogen gas. The calculations employ 
pseudopotentials and B-spline basis sets to expand the single particle orbitals used to construct the 
trial wavefunctions. Extensive tests on system size, time step, and other sources of errors, performed 
on periodically repeated systems of up to 1050 atoms, show that all these errors together can be 
reduced to below 10 meV per formula unit. We find excellent agreement with the experiments for 
the equilibrium volume of both the Mg and the MgH2 crystals. The cohesive energy of the Mg 
crystal is found to be 1.51(1) eV, and agrees perfectly with the experimental value of 1.51 eV. The 
enthalpy of formation of MgH2 from Mg bulk and H2 gas is found to be 0.85±0.01 eV/formula unit, 
or 82±1 kj/mole, which is off the experimental one of 76.1±1 kj/mole only by 6 kj/mole. This shows 
that DMC can almost achieve chemical accuracy (1 kcal/mole) on this system. Density functional 
theory errors are shown to be much larger, and depend strongly on the functional employed. 

I. INTRODUCTION 

The energetics of metal hydrides has recently become an issue of large scientific and technological interest, mainly 
because of the revived interest in these materials as potential hydrogen storage media^ . Magnesium hydride (MgH2) 
is a particularly interesting material, as it can store up to 7.6 % of hydrogen by weight, which is believed to be a large 
enough quantity for mobile applications, provided that all the hydrogen in the material can be made available when 
requested, of course. When heated above ~ 300 °CFMgH2 decomposes into Mg bulk and H2 gas, the reaction being 
endothermic with an enthalpy of decomposition of 76 kJ/mole'^. Conversely, MgH2 can be synthesised by combining 
Mg bulk (usually in form of a powder of micro- metre sized grains) and H2 gas. The charging process can take many 
hours, because of a large energy barrier to dissociate the H2 molecule on the surface of magnesium^. As it stands, 
MgH2 is not considered to be useful for hydrogen storage purposes, because of the high decomposition temperature 
(ideal decomposition temperature should be in the range 20 — 100 °C), and the slow kinetics of hydrogen intake. A 
number of attempt s are being made to modify this material to improve its properties, including doping it with traces 
of transitio n metalj^ *^, which have been shown to be very effective at reducing the activation energy for hydrogen 
dissociation ISinEnin] and also somewhat reduce the decomposition temperature of the hydride-. 

A number of theoretical calculations h ave been performed on magnesium hydride and related systems (see for 
example and references therein; see also I 6 | 8 | 9 | i0 | ii | i3 | i4 | i5 | i 6l ^^ ^-^^ most recent ones based on the implementation of 
quantum mechanics known as density functional theory (DFT/^^. Although DFT can often be reliable at predicting 
trends in the energetics of materials, it can be sometime in error when used to obtain absolute energies. In particular, 
as we show below, when applied to the calculation of the enthalpy of formation of MgH2, the results are off by as 
much as 0.3 eV per formula unit, depending on the functional employed, and cohesive energies can be wrong by over 
0.5 eV. 

Quantum Monte Carlo (QMC) technique^^^ are believed to be one possible way to improve beyond density 
functional theory. Since they are many order of magnitudes more computationally demanding, the current database 
of properties of materials calculated with QMC is still rather small, however, the increase in computer power in 
the past few years is making now possible to perform increasingly more numerous calculations on real systems, and 
experience is being accumulated on the predictive power of this technique. 
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Here we have used QMC to calculate the structural properties of the Mg and MgH2 crystals, together with their 
cohesive energies and the enthalpy of formation of MgH2 from Mg bulk and H2 gas. We find excellent agreement 
with experiments for the structural properties of the two solids, as well as the cohesive energy of the Mg solid. The 
enthalpy of formation of MgH2 is slightly overestimated, but the error is of the order of 1 kcal/mole, showing that 
QMC on this system can almost achieve chemical accuracy. 

II. TECHNIQUES 

A. Density Functional theory calculations 

Density functional theory calculations have been performed with the VASP cod^^H, Th e inte ractions between the 
electrons and the ionic cores was described using the projector augmented method (PAW^^SlSil with the generalised 
gradient approximations known as PBE^'"', PW91^^ or the local density approximation (LDA). The Mg PAW potential 
has a frozen Ne core and an outermost cutoff radius for the valence orbitals of 1.06 A. The H PAW potential has 
a cutoff radius of 0.58 A. Single particle orbitals were expanded in plane waves with a plane-wave cutoff of 270 eV, 
and a cutoff of 1600 eV was used for the charge density. Such a large cutoff^ in the charge density (4 times larger 
than the typical one used by default) is necessary to obtain very accurate forces which are used to calculate the 
vibrational properties of the crystals. Calculations were performed by requiring a self-consistency convergence on the 
total energy of 10~* eV per simulation cell. With these prescriptions convergence on the forces was at worse equal 
to 0.2 me V/A, and one or two order of magnitudes smaller for most atoms in the simulation cell. Brillouin Zone 
integration was performed using k-point sampling, with 18x18x12 and 10x10x15 Monkhorst-Pack^^ grids on the Mg 
and MgH2 primitive cells respectively. With these densities of k-points the structural parameters are converged to 
better than 0.1 %, and the total energies to better than 1 mc V/primitivc cell. 

B. Quantum Monte Carlo calculations 

Quantum Monte Carlo techniques have being extensively described elsewher d^^ l ^" !, so here we only report the main 
technical details used in this work. Calculations have been performed using the CASINO codtP. Diffusion Monte 
Carlo calculations have been performed using trial wavefunctions of the Slater-Jatrow type: 

^'t(R) = D^D^e' , (1) 

where and are Slater determinants of up- and down-spin single-electron orbitals, and e"^ is the so called 
Jastrow factor, which is the exponential of a sum of one-body (electron-nucleus), two-body (electron-electron) and 
three body (electron-electron-nucleus) terms, which are parametrised functions of electron-nucleus, electron-electron 
and electron-electron- nucleus separations, and are designed to satisfy the cusp conditions. The parameters in the 
Jastrow factor are varied to minimise the variance of the local energy El{1V) = ^'|p^(R)J?^'t(R). 

Imaginary time evolution of the Schrodinger equation has been performed with the usual short time approximation, 
and the locality approximation^®. Time step errors have been carefully analysed later in the paper. Since the locality 
approximation introduces an uncontrollable error with respect to which the DMC energy is non-variational, we also 
tested the scheme of CasulaP^, which treats the non local part of the pseudopotential in a consistent variational 
scheme. We found that the zero time step extrapolation of the energies in the two scheme differed very little, which 
suggests that the errors in either case is rather smalP^. However, we also found that the time step error is much 
smaller in the locality approximation in this particular case (this may not be true in general for other systems), and 
therefore we decided to use the locality approximation throughout the work which allowed us to work with a larger 
time step. 

We used Dirac-Fock pseudopotentials (PP) for Mg and lP3. The Mg PP has a frozen Ne core and a core radius 
of 1.43 A, the H PP has a core radius of 0.26 A. The single particle orbitals have been obtained by DFT plane-wave 
(PW) calculations using the LDA and a PW cutoff of 3400 eV, using the pwscf packagepS'. Such a large PW cutoff 
is due to the very small H PP core radius, and was found to be necessary to reduce the variance of the local energy 
as much as possible. We then exploited the approximate equivalence between PW and B-spline^^ to expand the 
single particle orbitals in a basis of B-splinc, as described in Ref.^'', using the natural B-spline grid spacing given by 
a — Tr/Gmax, where Gmax is the length of the largest vector employed in the PW calculations. 

We used a diffusion Monte Carlo time step of 0.05 a.u., which was found to result in errors of about 2 meV/f.u. (see 
below). With this time step the acceptance ratios were 99.2 and 99.7 % for the MgH2 and Mg crystals respectively. 
Total energies in the solids were obtained by correcting the raw DMC data with DFT-LDA calculations performed on 
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the same cell size but a fully converged Brillouin zone sampling, and then extrapolating these corrected DMC data to 
infinite size (see below) . The DMC calculations were performed using the Ewald interaction to model electron-electron 
interactions. The number of walkers in the DMC simulations varied with the size of the systems, and was never less 
than 1280. 



III. RESULTS 

Mg bulk has the hexagonal close packed structure, which is specified by a lattice parameter a and the ratio c/a 
of the vertical axis to one of the horizontal ones. The primitive cell contains two atoms, one at the origin and the 
other at (1/3,2/3,0.5) in lattice vectors units. The MgH2 solid has a tetragonal structure of rutile type (see Fig. [ij, 
specified by a lattice parameter a and the c/a ratio. The primitive cell has two Mg atoms, one at the origin and the 
other in the centre of the cell at (1/2,1/2,1/2) plus four hydrogen atoms at (±x,±x,0) and (l/2±x,l/2=Fx,l/2). The 
exact values of c/a and x depend on pressure, and at ambient conditions are found to be c/a = 0.6687 and x=0.304^. 

A. Zero point energies and high temperature vibrational effects 

In order to compare the calculated structural parameters and cohesive energies with the experimental ones we 
need to study the vibrational properties of the crystals. This is because the experimental parameters are usually 
determined at ambient conditions, and room temperature thermal expansion for the Mg and MgH2 solids is likely to 
be significant. 

We studied these vibrational properties within the quasi-harmonic approximation, which far from the melting 
temperatures provides accurate enough results for the thermal expansion of solids. This is certainly the case for the 
Mg and MgH2 solids at room temperature. 

Phonons have been calculated using the phon cod^^, which implements the small displacement methocP^I^ to 
obtain the force constant matrix in crystals. The methods exploits the linearity relation between the displacement of 
the atoms from their equilibrium positions and the forces induced on all the atoms in the crystal, which holds in the 
harmonic approximation for small enough displacements. The method is applied by constructing a supercell which is a 
multiple of the primitive cell in the three spacial directions, then the atoms in the primitive cell are displaced by small 
amounts along three linearly independent directions and the forces induced on all the atoms in the supercell are used 
to construct the force constant matrix. Symmetries can usually be used to reduce the total number of displacements 
needed, and also to symmetrise the force constant matrix:^^'. For bulk Mg, which has the hexagonal closed packed 
crystal structure, only two displacements are needed, one in the basal plane and one orthogonal to it (in fact, one 
single off symmetry displacements would be sufficient, although this would break the symmetry of the supercell and 
require a larger number of k-points in the DFT calculation of forces) . MgH2 has the tetragonal structure of rutile 
Ti02, with two Mg and four H atoms in the primitive cell, and the total number of displacements needed in this case 
is 4 (one could reduce the total number of displacements to 2 by sacrificing symmetries). If the supercell is large 
enough so that the forces on the atoms sitting near the edges are small, then the calculated force constant matrix 
becomes a good approximation of the exact one. Magnesium bulk is a metal, and convergence of the force constant 
matrix with the size of the supercell is readily achieved: we found that with cell containing 36 atoms (3x3x2) the 
ZPE is converged to within 0.1 meV/atom (tested using supercells containing up to 150 atoms). However, MgH2 is 
an insulator, and long range Coulomb interactions make convergence slower. Nevertheless, we found that already by 
using a cell containing 72 atoms (2x2x3 supercell) the ZPE can be calculated with an accuracy of 0.5 meV/fu (tests 
used supercells containing up to 576 atoms). All calculations were performed with DFT-PBE. 

Phonons calculated with the direct method described above may suffer from inaccuracies due to the size of the 
displacements and/or numerical noise in the calculated forces. To reduce the latter, one would like to maximise the 
size of the displacements, but too large displacements would cause departure from the harmonic regime. A compromise 
between these two opposite requirements then needs to be found, and this is usually achieved with displacement sizes 
of the order of a fraction of a percent of the inter-atomic distances. In order to test the size of the displacements we 
repeated the calculations using displacements of 0.067 A, 0.04 A, 0.02 A, and 0.01 A, and we found that even with 
the largest displacement the ZPE energy is converged to less than 0.2 meV/atom in Mg and 1 meV/fu in MgH2. We 
then decided to use displacements of 0.04 A. 

The fundamental vibrational frequency of the H2 molecule has been obtained by calculating the total energy of the 
H2 molecule in a large cubic box of size 13.5 A for 5 different values of the H-H distance, ranging from i?o-0.0135 Ato 
i?o -I- 0.0135 A, where Rq — 0.75 A is the calculated equilibrium distance with DFT-PBE. The 5 energies have been 
fitted to a parabola, providing a force constant of 33.35 eW / K? which corresponds to a stretching vibrational frequency 
of 127 THz (only shghtly lower that the experimental value of 131.8 THz^a), giving a ZPE of 0.263 eV. 
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TABLE I: Bulk properties (Volume/fu Vo in A^, and bulk modulus ko in GPa) and cohesive energies (Ecoh, in eV) of Mg 
and MgH2. Calculated properties are reported at zero temperature with and without zero point energies (ZPE) and at the 
temperatures at which the experimental data have been taken. Also reported is the binding energy of the H2 molecule 
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Vo, ko 


Mg 


LDA 
PBE 
PW91 

FT? .-, 1 

[Exp.J 
DMC 


21.59, 40.6 
22.86, 36.5 
22.86, 36.4 

22.96 ±0.05, 35.5 ± 1.2 


21.80, 39.3 
23.08, 35.9 
23.10, 35.2 

23.19 ±0.05, 34.4 ± 1.4 


-1.74 
-1.47 
-1.45 

F 1 C 1 1 

[-1.51] 
-1.51 ±0.01 


22.14, 36.4 
23.47, 34.0 
23.50, 32.6 

Foo O/IC oi^ 1 ci d 

[23.24 , 35.8 ± 3.0 
23.61 ±0.04, 31.2 ±2.4 


MgHs 


LDA 
PBE 
PW91 

FT? . , 1 

[Exp.J 
DMC 


29.36, 55.5 
30.84, 51.1 
30.72, 51.5 

29.48 ±0.03, 58.6 ±3.6 


30.32, 51.4 
31.92, 45.8 
31.79, 46.4 

30.53 ±0.05, 42.0 ± 1.5 


-7.16 
-6.17 
-6.27 
[—6.78 ± 0.01 ] 
-6.84 ±0.01 


30.36, 49.9 
32.03, 43.5 
31.89, 43.9 
[30.49 , - J 
30.58 ±0.06, 39.5 ± 1.7 


H2 


LDA 
PBE 
PW91 
[Exp] 
DMC 






-4.59 
-4.23 
-4.25 

[-4.48"] 
-4.484 ±0.002 





'^T= 298 K for Mg, T= 260 K for MgH2. 
''Ref. mi 

''Ret \M 
"Ret [1 
■''Ref. [35] 
sRef. [39] 

B. Density functional theory results 

Initially, we performed DFT calculations on the crystals with PBE, PW91 and LDA. Energy versus volume curves 
were fitted to a Birch-Murnaghan equation of state^, which provided equilibrium volumes and bulk moduli. In 
the range of volumes considered, c/a's do not change very much from their zero pressure values, and the structural 
parameters are essentially unchanged if c/a is kept fixed. Therefore, for simplicity we decided to fix c/a to their 
calculated zero pressure values of 1.621 and 0.6682 for Mg and MgH2 respectively. The MgH2 crystal has an additional 
degree of freedom, which defines the position of the H atoms in the lattice. This has also been optimised by fully 
relaxing the crystal at each different volume. These relaxations are essential in the calculation of phonons, because 
if the crystal is not in its ground state imaginary phonon frequencies appear. However, as far as the energy is 
concerned, the differences from calculations in which the H positions are kept at their zero pressure equilibrium values 
are undetectable. 

In Table |l] we report the structural parameters of Mg and MgH2 calculated with the three density functionals, and 
we report the results both at zero temperature (with and without ZPE) and at room temperature. Both Mg and 
MgH2 are fairly soft materials, with bulk moduli of the order of 40 and 50 GPa. Room temperature thermal pressure 
are about 1 and 1.8 GPa for Mg and MgH2 respectively, which means volume thermal expansion is about 2% and 3.5% 
for the two solids. This is significant, and cannot be ignored in a fair comparison with the experimental data. We also 
report in the same table the cohesive energies of the two solids. The experimental cohesive energy of MgH2 can be 
estimated by combining the cohesive energy of the Mg crystal (1.51 eV/atom), the dissociation energy of the hydrogen 
molecule (4.48 eV/molecule) and the enthalpy of formation of MgH2 from Mg and H2, whose value extrapolated at 
zero temperature is 0.79±0.01 eV/fiP, which therefore give a result of 6.78±0.01 eV/fu. By comparing the calculated 
cohesive energies with the experimental ones it is clear that the three functionals provide quite scattered results, with 
the LDA doing better on MgH2 and PBE doing better on Mg. It is also apparent that errors can be significant, of up 
0.6 eV for PBE. This error is well over 10 times a kcal/mole, which is the typical quantity cited as chemical accuracy. 
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C. Diffusion Monte Carlo results 

1. Time step tests 

The dependence of the DMC energy on time step in the MgH2 crystal was studied by repeating simulations with a 
2x2x3 supercell (72 atoms) at time steps ranging from 0.005 to 0.15 a.u.. Calculations were performed at the volume 
of 30.835 AVfu, and using the A point (0.5,0.5,0.5) which is at one corner of the Brillouin zone. For the Mg crystal 
we used a 3x3x2 supercell (36 atoms), a volume of 22.785 A'^/atom and the H point (0.5,0.5,0.5), also at one corner 
of the Brillouin zone. 

Results of total energy/fu for MgH2 and total energy/atom for Mg are displayed in Fig. [2] from which it is evident 
that using a time step of 0.05 a.u. time step errors are well below 5 meV/fu. In Fig. |2]we also display the results 
obtained with the scheme proposed by CasulaP^]^ and we observe that for short enough time steps the two sets of 
energies are very close, and extrapolate to roughly the same value in the limit of zero time step (to less than 5 
meV/fu). As mentioned earlier, this suggests that the error introduced with either scheme is very small. However, 
the locality approximation results in a much weaker dependence of the DMC energy on time step and this is what we 
used because it allowed us to work with much larger time steps. We note that for the Mg crystal the time step error 
is much smaller, which in principle would allow us to work with larger time steps, however, for consistency, we used 
the same time step of 0.05 a.u. also for the Mg crystal. 

To calculate the total energies of the Mg atom and the H2 molecule we used trial wavefunctions obtained from 
plane wave calculations in which the Mg atom or the H2 molecule was placed at the centre of a large cubic box with a 
side of 13.5 A. The DMC calculations were then performed using B-splines and no periodically boundary conditions. 
We display in Fig. [3] the DMC energies as function of time step, from which we can obtain very accurate zero time 
step values. In the case of Mg we also performed one calculation with the scheme of CasulaP^, which gave essentially 
the same energy. For the H2 molecule we display the binding energy calculated at the equilibrium distance of 0.75 A, 
obtained by subtracting from the energy of the molecule twice the energy of the H atom, which is calculated to be 
13.60635(5) eV. Both the energies of the H atom and the H molecule are in excellent agreement with the experimental 
data. 



2. The Mg crystal 

In the Mg crystal we studied the dependence of the DMC energy on the size of the simulation cell by repeating 
the calculations with 4x4x3, 5x5x3, 6x6x4, 8x8x5 and 9x9x6 supercells, containing 96, 150, 288, 640 and 972 atoms 
respectively. Results are displayed in Fig.|4] where we show the total energies/atom Ejq as function of 1/A^, with N the 
number of atoms in the simulation cell. On the same graph we also show the energies i?^ = Etq + \E^^'^ — i?^^^] , 
where E^^^ are the DFT energies calculated with fully converged k-point sampling, and E^^'^ are the DFT energies 
calculated with k-point samplings corresponding to the A^-atom cells used in the DMC calculations. It is clear that 
the raw DMC energies E^ are quite scattered and somewhat difficult to extrapolate to infinite size. This is due to 
the metallic nature of Mg. However, the DFT corrected energies E"^ are much better behaved, with data fitting quite 
well onto a straight line, which makes it possible to extrapolate to infinite size. In particular, we note that with no 
loss of accuracy we can also use only the calculations with the 4x4x3, 5x5x3 and 6x6x4 supercells to extrapolate to 
essentially the same infinite size value. 

The calculations with these three supercell sizes were then repeated at 8 different volumes, between 21 and 
25 A^/atom. At each volume the DFT corrected DMC results were extrapolated to infinite size and the results 
were fitted to a Birch-Murnaghan equation of state to obtain the structural parameters. We performed the fit by 
weighing each energy point Ei point with l/af , where Ui is the standard error on Ei. We report in Table |l] the results 
obtained both at zero temperature (with and without zero point energy) and at room temperature. The latter are also 
shown in Fig. [5] The room temperature corrected DMC results slightly overestimate the equilibrium volume, and also 
underestimate the bulk modulus, but the calculated cohesive energy is in perfect agreement with the experimental 
data. 



3. The MgH2 crystal 

For the MgH2 crystal size effects were studied using 2x2x3, 3x3x4, 4x4x6 and 5x5x7 supercells, containing 72, 216, 
576 and 1050 atoms respectively. These tests were performed at the volume of 30 A^^/fu. The results for the four sizes 
studied are displayed in Fig. [6] where we show total energies/fu E'at as function of 1/A, as well as the DFT corrected 
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energies i?^. In this case the DFT corrections are much smaller, which is not surprising because of the large band 
gap in MgH2. A small difference between the two sets of data can be observed for the smallest sizes, but it is clear 
that they both fit very well onto straight lines, which allows us to easily extrapolate the results to infinite size. In 
fact, in this case the extrapolated results for the two sets only differ by 5 meV/atom. 

The calculations were repeated at 7 different volumes between 28 and 32.5A'^/fu, and the DFT corrected DMC 
results were then fitted with a Birch-Murnaghan equation of state to obtain the structural parameters. Also in this 
case we used the inverse of the variances to weigh each point in the fit. We report in Table [T] the results obtained 
both at zero temperature (with and without zero point energy) and at T = 260 K, which is the temperature at which 
the experimental data are reportecP^. The high temperature results are also shown in Fig. [7] It is clear that once 
thermal effects are added onto the calculations the agreement with the experimental equilibrium volume is in almost 
perfect agreement. The cohesive energy is slightly overestimated, but the error is only 0.06 eV, i.e. of the order of 
chemical accuracy. 

4- Enthalpy of formation of MgH2 

We can now calculate the enthalpy of formation of Mgll2 from Mg bulk and H2 in the gas phase by adding the 
cohesive energies of the MgH2 and Mg crystals to the binding energy of the II2 molecule. We obtain enthalpy of 
formations of 0.82, 0.47 and 0.57 eV/fu with LDA, PBE and PW91 respectively, and with DMC we obtain the value 
0.85± 0.01 eV/fu. The LDA value is very accurate, but this is the result of large cancellations of errors in the cohesive 
energies of the crystals and the binding energy of the II2 molecule. The DMC result is only 0.06 eV higher than 
the experimental value of 0.79 ± 0.01 cV/fu, however in this case both the cohesive energies of the crystals and the 
binding energy of the H2 molecule are very accurate. 

IV. CONCLUSIONS 

We pointed out in this work the difficulty of using density functional theory to calculate the enthalpy of formation 
of MgH2 with high accuracy. We studied the effect of three different exchange-correlation functionals, PW91, PBE 
and LDA, and found that although the GGA ones appear to work better on the Mg solid, the LDA gives better results 
on the MgH2 solid. It turns out, therefore, that is difficult to get a good DFT value for the enthalpy of formation of 
MgH2: the two GGA functionals give an enthalpy of formation in error of more than 0.2 and 0.3 eV/fu respectively. 
The LDA is the functional that does best, but for the wrong reason, because the cohesive energies of the crystals and 
the binding energies of the molecule are wrong by up to 0.4 eV/fu, and the enthalpy of formation is accurate only 
because of large cancellation of errors. 

Diffusion Monte Carlo appears to deliver much better accuracy in general. We have shown that the DMC equilibrium 
volumes of MgH2 agrees perfectly with the experimental one, once high temperature thermal expansion is included 
in the calculations, and the equilibrium volume of Mg is only slightly overestimated. The cohesive energy of Mg is 
also predicted in perfect agreement with the experimental datum, and so is the binding energy of the H2 molecule. 
A small error is present in the cohesive energy of the MgH2 crystal, which determines the small inaccuracy in the 
enthalpy of formation, for which we find a DMC value of 0.85 ± 0.01 eV/fu. However, this is only 0.06 higher than 
the accepted experimental one of 0.79 ± 0.01 eV/fu, or 76.1 ± 1 kJ/mole. This result is not very far from the LDA 
value, but with the important difference that now all three terms that enter the enthalpy of formation are calculated 
accurately, and we don't rely on fortunate cancellation of errors. 

Although the DMC error is slightly larger than 1 kcal/mole, and therefore we cannot claim chemical accuracy, we 
are not far from it, and therefore we argue that quantum Monte Carlo techniques have useful predictive power in the 
search of metal hydrides with workable decomposition temperatures. 

Acknowledgments 

The computations were performed on the HPCx service, using allocations of time from NERC, and on HECToR. 
Calculations have also been performed on the LCN cluster at University College London. This work was conducted 



7 



as part of a EURYI scheme award as provided by EPSRC (see www.esf.org/euryi). 



2 



Electronic address: d.alfe @ucl.ac. uk ' 

L. Schlapbach and A. Ziittel, Nature (London) 414, 353 (2001). 

B. Bogdanovic, K. Bohmhammel, B. Christ, A. Reiser, K. Schlichte, R. Vehlen and U. Wolf, J. Alloys Comp. 282, 84 (1999). 
^ M. Yamaguchi and E. Akiba, in Material Science and Technology, vol. 3B, edited by R. W. Calm, P. Haasen and E. J. 

Kramer (New York: VCH, 1994), p. 333. 

P. T. Sprunger and E. W. Plummer, Chem. Phys. Lett. 187, 559 (1991). 
^ G. Liang, J. Huot, S. Boily A. Van Neste and R. Schulz, J. Alloys Compd. 292, 247 (1999). 
^ C. X. Shang, M. Bououdina, Y. Song and Z. X. Guo, Int. J. Hydrogen Energy 29, 73 (2004). 
^ N. Hanada, T. Ichigawa and H. Fujii, J. Phys. Chem. B 109, 7188 (2005). 
^ A. J. Du, S. C. Smith, X. D. Yao, and G. Q. Lu, J. Phys. Chem. B 109, 18037 (2005). 
® A. J. Du, S. C. Smith, X. D. Yao, and G. Q. Lu, J. Phys. Chem. B 110, 21747 (2006). 
A. J. Du, S. C. Smith, X. D. Yao, and G. Q. Lu, J. Am. Chem. Soc. 129, 10201 (2007). 
M. Pozzo, D. Alfe, A. Amieiro, S. French and A. Pratt, unpublished. 
12 R. Yu and P. K. Lam, Phys. Rev. B 37, 8730 (1988). 

1^ P. Vajeeston, P. Ravindran, A. Kjekshus, and H. Fjellvag, Phys. Rev. Lett. 89, 175506 (2002). 
" Y. Song, Z. X. Guo and R. Yang, Phys. Rev. B 69, 094205 (2004). 

1^ M. J. van Setten, G. A. de Wijs, V. A. Popa, and G. Brocks, Phys. Rev. B 72, 073107 (2005). 

T. Vegge, L. S. Hedegaard- Jensen, J. Bonde, T. R. Munter and J. K. N0rskov, J. Alloys Compd. 386, 1 (2005). 
^'^ P. Hohenberg, W. Kohn, Phys. Rev. B 136, 864 (1964). 
1* W. Kohn and L.J. Sham. Phys. Rev. A 140, 1133 (1965). 

1® W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001). 

C. J. Umrigar, M. P. Nightingale, K. J. Runge, J. Chem. Phys., 99, 2865 (1993). 

R.J. Needs, M.D. Towler, N.D. Drummond and P. Lopez Rios, CASINO version 2.1 User Manual, University of Cambridge, 
Cambridge (2007). 

22 G. Kresse and J. FurthmuUer, Comput. Mater. Sci. 6 (1996) 15; Phys. Rev. B 54, 11169 (1996). 
P. E. Blochl, Phys. Rev. B, 50, 17953 (1994). 

G. Kresse and D. Joubert, Phys. Rev. B, 59, 1758 (1999). 
J. P. Perdew, K. Burke and M. Ernzerhof. Phys. Rev. Lett. 77, 3865 (1996). 
J. P. Perdew and Y. Wang. Phys. Rev. B 45, 13244 (1992). 

H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976). 
L. Mitas, E. L. Shhley D. M. Ceperley J. Chem. Phys. 95, 3467 (1991). 
M. Casula, Phys. Rev. B 74, 161102 (2006). 

Of course, it is possible that the Casula scheme induces a large positive error in the DMC energy, and that by chance 
the error induced by the locality approximation is also positive and of the same size. However, since the two scheme are 
completely different, we believe that this is unlikely. 

J. R. Trail and R. J. Needs, J. Chem. Phys. 122, 014112 (2005); ibid J. Chem. Phys. 122, 174109 (2005), see also 
www.tcm.phy.cam.ac.uk/~mdt26/casino2_pseudopotentials.html. 



11 



20 



32 
33 
34 



S. Baroni, A. Dal Corso, S. de Gironcoli, and P. Giannozzi, http://www.pwscf.org 
E. Hernandez, M. J. Gillan and C. M. Goringe, Phys. Rev. B 55, 13485 (1997). 
D. Alfe and M. J. Gillan, Phys. Rev B, 70, 161101(R) (2004). 
M. Bortz, B. Bertheville, G. Bottger, K. Yvon, J. A lloys Comp d. 287, L4 (1999) 

D. Alfe, (1998). Program available at http://chianti.geol.ucl.ac.uk/~darioi D. Alfe, G. D. Price, M. J. Gillan, Phys. Rev. B, 
64, 04512316 (2001). 



G. Kresse, J. Furthmiiller and J. Hafner, Europhys. Lett. 32, 729 (1995). 
D. Alfe, G. D. Price, M. J. Gillan, Phys. Rev. B, 64, 045123 (2001). 

B. H. Bransden and C. J. Joachin, Physics of Atoms and Molecules (Wiley, New York, 1983). 
F. Bhch, Phys. Rev. 71, 809 (1947). 
*i C. Kittel, Introduction to Solid State Physics, 7th ed. (Wiley, New York, 1996). 
*2 http://www.webelements.com 

D. Errandonea, Y. Meng, D. Hausermann, and T. Uchida, J. Phys.: Condens. Matter 15, 1277 (2003). 



8 




Mg bulk MgHj bulk 




dt (a.u.) dt (a.u.) 



FIG. 2: Diffusion Monte Carlo energies for Mg bulk (left panel) and MgH2 bulk (right panel) as function of time step. Dots 
and squares correspond to calculations performed with the locality approximation and with the scheme proposed by Casul^^ 
respectively. 
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Mg atom molecule 




dt (a.u.) dt (a.u.) 



FIG. 3: Dots: diffusion Monte Carlo total energy for the Mg atom (left panel) and binding energy of the H2 molecule (right 
panel) as function of time step. Calculations have been performed with the locality approximation. Square: calculation 
performed with the scheme proposed by Casuld^. 
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FIG. 4: Diffusion Monte Carlo total energy for the Mg crystal as function of 1/N, where N is the number of particles in the 
simulation cell. Stars and squares correspond to raw and DFT corrected (see text) results, solid line is a linear least square fit 
to the DFT corrected results. 
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FIG. 5: Diffusion Monte Carlo free energies at 298 K for the Mg crystal as function of volume V. Dots correspond to DMC 
calculations extrapolated to infinite size, and include vibrational free energies calculated with DFT-PBE. Solid line is a least 
squares fit to a Birch-Murnaghan equation of state. 
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FIG. 6: Diffusion Monte Carlo total energy for the MgH2 crystal as function of where N is the number of particles in the 
simulation cell. Stars and squares correspond to raw and DFT corrected (see text) results, solid line is a linear least square fit 
to the DFT corrected results. 
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V (A'/f.u.) 

FIG. 7: Diffusion Monte Carlo free energies at 260 K for the MgH2 crystal as function of volume V. Dots correspond to DMC 
calculations extrapolated to infinite size, and include vibrational free energies calculated with DFT-PBE. Solid line is a least 
squares fit to a Birch-Murnaghan equation of state. 



